Inter-annual variations of vegetation dynamics to climate change in Ordos, Inner Mongolia, China

To reveal the characteristics of climate change and the controlling factors for vegetation dynamics in the Ordos, Inner Mongolia, China, 34 years (1982–2015) of regional climate variables and vegetation dynamics were investigated. The results show that: Annual mean air temperature (TMP) significantly increased with a linear slope of 0.473°C/10yr. Annual precipitation (PRE) had a non-significant positive trend nearly 5 times lower than the trend of potential evapotranspiration (PET). The average Normalized Difference Vegetation Index (NDVI) computed for the region was found to show a significant positive trend (6.131×10−4/yr). However, all climate variables displayed non-significant correlations with NDVI at annual scale. The reduction of desert and the increase of grassland over the past decades were accountable for the increased NDVI. Principal components analysis revealed that the regional climate change can be characterized as changes in temperature, humidity and the availability of radiant energy. Based on principal components regression coefficients, NDVI was mostly sensitive to humidity component, followed by growing season warmth (WMI). Spatially, 93.1% of the pixels displayed positive trend and 61.8% of the pixels displayed significant change over the past decades. Both principal regression analysis and partial correlation analysis revealed that NDVI in eastern part of Ordos was sensitive to TMP, whereas, NDVI in southern and western areas of Ordos displayed the high sensitivity to combined effects of PRE and cloud coverage (CLD). Partial correlation analyses also revealed that TMX was a surrogate for aridity, TMN was a representative of humidity, and temperature variations below the threshold of 5°C (CDI) were less important than WMI. We conclude that regional climate change can be characterized by warming and increased aridity. The significant positive trend of regional NDVI and the non-significant correlations between NDVI and climate variables at annual scale suggests the hidden role of the human activities.


1.Introduction
Vegetation is a critical component of the terrestrial biosphere and its growth activity plays an important role in the global carbon and hydrological cycles [1]. Studies have shown that there are clear relationships between climate data and Global Inventory Modeling and Mapping Studies (GIMMS) satellite-derived normalized difference vegetation index (NDVI) vegetation dynamics [2]. Additionally, feedbacks have been revealed between vegetation activity and the carbon and hydrological cycles [3][4][5][6][7].
Most of these studies concentrated exclusively on empirical data such as observed temperature and precipitation. Generally, these studies do not consider derived climate variables such as evapotranspiration [8], cloud cover [9], and warmth/coldness indices [10]. In addition, the effects of asymmetrical climate change, which is characterized as a higher increase of minimum temperatures and a smaller increase in maximum temperatures, on vegetation growth have not been sufficiently investigated [11,12].
Traditionally, studies that have identified significant climate-vegetation relationship with linear correlation analytical approach [13]. Limited consideration has been given to the co-linearity among climate variables, which may results in pseudo driving factors for vegetation dynamics [14,15].
Moreover, it is challenging to differentiate between the climate and human causations on vegetation dynamics. Urbanization, cropland expansion, demographic changes, and livestock grazing activity have been found to impact climate-vegetation dynamics, especially within the forest-steppe-desert transition zone where different trends and driving factors coexist [6,16,17]. Traditional linear correlation approaches provide an incomplete representation of variable-vegetation analysis, which may lead to biased conclusions on the scope of climate impacts on vegetation and thereby influence regional climate-vegetation modelling scenarios. For example, the expansion of cropland in north-eastern Inner Mongolia misled our understanding about the climate-vegetation relationship [10]. A climate-vegetation study in southern Africa underestimated the severity of land degradation using the linear regression method [18].
The Ordos Plateau of Inner Mongolia, China, is a typical forest-steppe-desert transition zone. The ecosystem diversity and vegetation dynamics in this region are highly sensitive to both climate change and anthropogenic influences [19]. These climate-vegetation relationships are complicated with asymmetric climate change and human activities that take place over the past three decades. These relationships are further complicated by numerous national and local ecological restoration projects that had been implemented by the government since the year 2000 [20,21]. Therefore, investigations on the climatic variations and human activities on regional vegetation changes are critical to characterize the vegetation growth responses to climate change. As underlying causes are identified in a robust manner, further progress can be made on the restoration and reconstruction policies of the Ordos Plateau ecosystems, which could be expanded into the vicinity regions with similar climate settings.
In this study, climate datasets were compared to vegetation dynamics in the Ordos Plateau over a 34-year period (1982-2015) to evaluate the regional climate-vegetation relationship. This study aims to i) characterize the temporal patterns during the regional climate change, ii) evaluate spatial patterns of vegetation-climate relationships, and iii) identify the controlling factors that drive vegetation dynamics in the Ordos ecosystems.

Study site
The Ordos region is located in the southern part of the Inner Mongolia Autonomous Region (37˚41 0 -40˚51 0 N; 106˚42 0 -111˚31 0 E) and occupies approximately 87,400 km 2 [22]. The area is bounded on the north by the Hobq Desert and the south by the Mu Us Sandy Land [22,23]. The major ecosystems are grassland, sandy land, cropland, meadows and forested land [22,24].
The study site experiences a temperate continental climate with high evaporative demand during the growing season, which typically lasts from May to October. The mean annual temperature is 5.3-8.7˚C with a high of 23.5˚C during the hottest month (July) and a low of -11.5˚C during the coldest month (January). The total annual precipitation is 150-450 mm with more than 70% of the total annual precipitation occurring from July to September. Generally, air temperature declines with elevation, whereas, precipitation decreases from east to west across the region [13]. Geologically, the region is sub-divided into the northern Yellow River alluvial plain, the eastern hilly region, the western arid plateau and the southern sandy plain; these geological features influence the vegetation types and the strength of aridity.

Data collection
To capture the temporal dynamics of vegetation change in response to climate change, 34 years (1982-2015) of Normalized Difference Vegetation Index (NDVI) of Ordos were obtained from the Global Inventory Monitoring and Modeling System (GIMMS) project (https://ecocast.arc.nasa.gov/data/pub/gimms/, accessed at 10/Oct/2018). To account for the issues of orbital drift, sensor degradation and radiation effect of volcanic eruption, GIMMS NDVI3g data have been normalized to produce a non-stationary 1981-2012 AVHRR NDVI3g time series [25]. The spatial resolution of the data was at 1/12˚and pixels with low quality data (flag value = 4-7) were removed from the analysis [26].
To reveal the characteristics of regional climate change, climate datasets (Table 1) were obtained from the Climatic Research Unit (CRU: University of East Anglia Climatic Research Unit, UK) data archives [27]. These CRU datasets contain homogenized monthly time series of precipitation, daily maximum and minimum temperatures, cloud cover, and other variables, all of which were gridded to 0.5x0.5 degree resolution. To match the spatial-temporal resolutions of GIMMS ndvi3g and facilitate the calculation at pixel scale, the climate datasets (0.5˚×0.5˚resolution) were disaggregated into 1/12˚by bilinear interpolation of the neighboring pixels.
As warmth index (WMI) and coldness index (CDI) were effective indicotors for vegetation growth and distribution limits [10,15,28], we followed You et al. [10] to established annual anomalies of WMI and CDI at each pixel by counting the annual sum of positive and negative differences between monthly means and 5˚C [29].
To obtain further information on vegetation trends, MODIS land cover product (MCD12Q1) for 2001 and 2015 were collected from https://e4ftl01.cr.usgs.gov/MOTA/ MCD12Q1.006/. The definition of land use types followed the International Geosphere-

PLOS ONE
Biosphere Program (IGBP) land cover classification system, and the accuracy of the land cover dataset has been validated by field observations [30].

Analytical methods
To reveal the long-term trends of climate variables and NDVI, both non-parametric trend tests and linear trend tests were conducted [31,32]. To reduce the influence of autocorrelation on trend detection, the datasets were pre-whitened (MK-TFPW) prior to applying the Mann-Kendall trend test [33]. To reduce the influence of abnormal anomalies and outliers on the results of the trend analyses, the non-parametric Sen [34] slope was also obtained by calculating the median slope between all pairwise combinations of points in time series of climate variables and NDVI. To reveal any significant change points within time series of climate variables and NDVI, the Mann-Whitney-Pettit test was conducted [35]. Detailed time series analyses and correlation analyses were presented by You et al. [10].
To remove the impact of co-linearity between the climate variables, principal components regression (PCR) was used to identify the relative importance of each variable driving interannual variations of NDVI [9,10,15]. In each pixel, we first extract 3 principal components from the climate dataset, representing temperature, humidity and global radiation components. Then, we established multi-linear regression model by using standardized NDVI as dependent and using the loading scores of the 3 components as independents. Then, we multiplied the loadings of each variable by the aforementioned multi-linear regression model coefficients and summed these scores. This enabled us to estimate the relative importance of each variable in driving the interannual variations of NDVI. For the annual means of regional climate variables and NDVI values, we extract 6 principal components representing more than 99% of variance of the climate dataset.
The partial correlation was conducted to reveal the correlation between X and Y on condition of the effect due to the third variable Z was eliminated (Eq 1). Suppose we are given three variables X 1 , X 2 and X 3 . Let r 12 denote the correlation coefficient between the variables X 1 and X 2 . Then the partial correlation coefficient between X 1 and X 2 by controlling for the effect of X 3 is denoted by the symbol r 12,3 . It is given by the following formula: r 12;3 ¼ r 12 À r 13 r 23 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ð1 À r 2 13 Þð1 À r 2 23 Þ p ð1Þ

Trend analyses
Regional climate change and vegetation dynamics for the Ordos region during the 1982-2015 study interval are displayed in Table 2 and Fig 1. Significant change points were detected during the year 1996-1997 in all of the temperature-related and cloud cover variables (TMP, TMX, TMN, WMI, CDI and CLD). Temperature-related variables displayed an increasing inter-annual trend prior to the 1996 change point (TMP, TMX, WMI and CDI). These variables then displayed a negative trend after 1996. At decadal scales, TMP showed an increasing slope (0.473˚C/10yr) (Fig 1A), but the increase was not significant. Prior to the change point, TMN significantly increased (0.567˚C/10yr; p<0.01) (Fig 1B), while TMX increased but was nonsignificant (0.473˚C/10yr) (Fig 1C). These two trend trajectories suggest that asymmetric warming had occurred and that minimum temperatures were rising faster than maximum temperatures. PRE displayed a nonsignificant positive trend (5.5 mm/10yr; nonsig.) (Fig 1F). Conversely, PET displayed a significant positive trend (29.4 mm/10yr; p<0.01) (Fig 1G), which suggests increased aridity took place in the region. As a consequence, VAP trends significantly decreased (-0.205 hPa/10yr; p<0.01) with a significant change point in the year 2003 ( Fig 1H). The regional NDVI displayed a significant positive trend (6.131×10 −3 /10yr; p<0.05) with a significant change point in the year 2006 ( Fig 1J)

Principal component analysis
For the inter-annual means of climate variables, principal component analysis shows the 1st principal component of the 9 climate variables was highly correlated with TMP, TMN, TMX, WMI, CDI and PET; all of which represent the temperature component ( Table 3). The 2nd principal component was highly correlated with PRE and VAP, which together represent the humidity component. The 3rd principal component was highly correlated with CLD, which representing the solar radiation component. Therefore, regional climate change can be summarized as the changes in temperature, humidity and the availability of radiant energy.
For the inter-annual means of NDVI, the correlation analyses showed that all climate variables displayed non-significant correlations with NDVI (Table 3). This lack of correlation suggests regional vegetation dynamics were not controlled by individual climate variables exclusively. The principal regression coefficients (PRC) show that NDVI was mostly sensitive to humidity component, followed by WMI. This sensitivity indicates humidity and growing season warmth were important control factors on the inter-annual variation of NDVI.

Spatial distribution of NDVI trend
The majority of the pixels (93.1%) displayed a positive trend in NDVI during 1982-2015 (Fig  2A). Additionally, most of the pixels (61.8%) displayed a significant change in the trend of NDVI over the past decades (Fig 2B). Principal regression analyses revealed that the change in vegetation in the eastern part of Ordos was sensitive to TMP. In the southern and western areas of Ordos, vegetation change was sensitive to the combined effects of PRE and CLD ( Fig  2C). Overall, the R 2 from principal regression analysis, residuals of climate-vegetation relationship, was low, suggesting the weak link between climate and vegetation and the hidden role of anthropogenic influences (Fig 2D).

Land use and land cover changes
The land use and land cover transition matrix from 2001 to 2015 displayed that the groundcover area of grasslands, shrubland, urban lands and croplands had increased from 2001 to 2015 (Table 4). Remarkable changes are the the increase of grassland and the decrease of desert land. Additionally, this analysis showed that 21.1% of the total desert area in 2001 was replaced by grassland in 2015.

Partial correlation analysis
Partial correlation analyses were conducted to reveal the comparative roles of individual climate variables (TMP vs. PRE, TMN vs. TMX, WMI vs. CDI) on the inter-annual variation of NDVI (3). NDVI of eastern and south-eastern parts of the region were positively correlated with TMP (controlling PRE) (Fig 3A), whereas, the rest of the region displayed positive correlation with PRE (controlling TMP, Fig 3B). Conversely, negative correlation with TMP (controlling PRE) were identified for the western part of the region. This negative correlation reveals the effects of drought had combined with increased TMP and PET that resulted in restricted vegetation growth. The analysis also reveals that elevation has an influence on which climate variables were important controls on NDVI. NDVI in low elevations in the eastern and southern parts of the region were negatively correlated with TMX (controlling TMN), but were positively correlated with TMN (controlling TMX) (Fig 3C and 3D). NDVI in the higher elevations in the central part of the region were positively correlated with TMX (controlling TMN), but negatively correlated with TMN (controlling TMX). Partial correlation between WMI and NDVI (Fig 3E, controlling CDI) displayed similar spatial pattern with Fig 3A, suggesting the temperature variations below the threshold of 5˚C is less important in controlling the vegetation growth. The

PLOS ONE
Joint effects of climate change and human activities on regional vegetation dynamics.
insignificant role of coldness on controlling the vegetation growth was confirmed by non-significant partial correlation coefficients between CDI and NDVI (controlling WMI, Fig 3F).

Characteristics of climate change
This study expanded on the traditional understanding of regional climate warming, it also revealed asymmetric climate warming characterized by higher trend in TMN than TMX. This disparity results in decreased diurnal temperature range (TMX-TMN) [12,36]. As the daily temperature minima is influenced by night time net radiation, increased TMN was likely caused by increased night time cloud coverage. This observation coincides well with the positive trend of CLD found in this study ( Table 2). The increase in CDI also coincided with winter

PLOS ONE
warming, which has been widely reported [37]. The thermal change points that took place in 1996 coupled with decreased CDI during post-1996 period coincide with the recovery of Siberian High intensity [38,39]. The significant change point of temperature-related variables may also be related to 1997/1998 El Niño event [40][41][42][43]. This study revealed a positive but nonsignificant trend of PRE which is coincided with previous studies within the Loess Plateau [13] and in Northwest China [17]. Increased PRE may have also driven increased actual ET through vegetation transpiration [21,44]. However, the benefits of positive PRE trend are negated by significant rising PET, which was largely thermically driven. Higher PET caused the regional climate to shift towards warming and drying trends. Consequently, the regional climate change could be summarized as warming and drying and influenced by larger scale climate dynamics such as the El Niño, Southern Oscillation and the Siberian High intensity.  [19,45] and sustainable development [24,46], resulting in the positive trends of vegetation growth and regional actual evapotranspiration [21]. As a result of high variance of NDVI in Ordos [47], annual scale climate variables displayed non-significant correlation with NDVI (Table 3), which suggests there is a weak link between climate change and vegetation dynamics over the past three decades [23]. Similarly, R 2 values of the principal regression analysis were low (Fig 2D), suggesting the limited contribution of inter-annual climate variability on vegetation dynamics. Previous studies also reported that the regional NDVI displayed opposite trend with the prediction using climate variables as model inputs [23,48]. Moreover, the vegetation resilience and the spatial distribution of vegetation types were accountable for the non-significant inter-annual climate-vegetation relationship (Table 5). For example, half of the plateau displayed positive correlation between precipitation and NDVI, whereas, weak correlation coefficients were mainly distributed in Hobq Desert and the Mu Us Sandy Land [22].
Besides, some edaphic and topographic factors were accountable for the spatial distribution of NDVI [53,54]. For example, the lithology of bedrock greatly affects vegetation cover and distribution in the Mu Us Sandy Land area, whereas, a high percentage farmlands and grasslands with large NDVI values are mainly distributed on low-permeability strata, such as the Quaternary Lake and alluvial deposits [49]. Moreover, it is also likely that vegetation dynamics may links to depletion of ground water [21], possibly leading to the rapid loss of regional lakes [55].
To further address the influence of human activity on climate-vegetation relationship, two typical pixels with negative NDVI slope were selected. Land use change from 1984 to 2015 at the 1st pixel (39˚41 0 50.6@N, 109˚56 0 20.45@E) displays the expansion of urban land and the 2nd pixels were related to lake extinctions (40˚09 0 44.46@N, 108˚27 0 28.04@E). In addition, the Global Human Influence Index Dataset of the Last of the Wild Project, Version 2, 2005 (LWP-2), was used as a surrogate for human activities (data not shown). This dataset derived from nine Table 5. Available studies on regional vegetation dynamics and its controlling factors.

Authors
Temporal Span

Data Sources Vegetation Trend Controlling Factors Human effects
Ma et al. [

PLOS ONE
global data layers covering human population pressure (population density), human land use and infrastructure (built-up areas, nighttime lights, land use/land cover), and human access (coastlines, roads, railroads, navigable rivers) [56] and had been applied to study the influence of human activity on species distribution, biodiversity protection and vegetation activities [15,57]. Preliminary investigation on the correlation between the the spatial distribution of human influence index and NDVI trends suggests the positive role of human influences. Further studies are needed to address the detailed information of the influence of human activity on vegetation growth and dynamics.